library(tidyverse)
library(cowplot)
********************************************************
Note: As of version 1.0.0, cowplot does not change the
default ggplot2 theme anymore. To recover the previous
behavior, execute:
theme_set(theme_cowplot())
********************************************************
library(lubridate)
Attaching package: ‘lubridate’
The following object is masked from ‘package:cowplot’:
stamp
The following objects are masked from ‘package:dplyr’:
intersect, setdiff, union
The following objects are masked from ‘package:base’:
date, intersect, setdiff, union
library(mgcv)
source("UVP_2017_library.R")
theme_set(theme_cowplot())
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
bes<- read_csv("dataOut/binned_EachSize.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
bds <- read_csv("dataOut/binned_DepthSummary.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
ues <- read_csv("dataOut/unbinned_EachSize.csv")
Parsed with column specification:
cols(
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m,
depth = [32mcol_double()[39m,
psd_gam = [32mcol_double()[39m,
vol = [32mcol_double()[39m,
sizeclass = [31mcol_character()[39m,
lb = [32mcol_double()[39m,
ub = [32mcol_double()[39m,
binsize = [32mcol_double()[39m,
TotalParticles = [32mcol_double()[39m,
nparticles = [32mcol_double()[39m,
n_nparticles = [32mcol_double()[39m,
biovolume = [32mcol_double()[39m,
speed = [32mcol_double()[39m,
flux = [32mcol_double()[39m,
flux_fit = [32mcol_double()[39m,
GamPredictTP = [32mcol_double()[39m
)
uds <- read_csv("dataOut/unbinned_DepthSummary.csv")
Parsed with column specification:
cols(
.default = col_double(),
project = [31mcol_character()[39m,
profile = [31mcol_character()[39m,
time = [34mcol_datetime(format = "")[39m
)
See spec(...) for full column specifications.
PhoticBase <- 160
OMZBase <- 850
PlotNParticles <- uds %>%
ggplot(aes(x = tot_nparticles, y = depth, col = profile)) +
facet_wrap(~project) +
geom_point(alpha = 0.3, shape = 1) +
scale_y_reverse() + scale_x_log10()
PlotNParticles
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(tot_nparticles~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG2 <- gam(tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG3 <- gam(tot_nparticles ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
#FSG4 <- gam(tot_nparticles~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4,
bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.002 0.106 84.96 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.204 1.366 3.346 0.0874 .
s(Day) 1.512 1.761 3.001 0.0386 *
s(Hour) 1.579 2.000 2.713 0.0389 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.183 Deviance explained = 24.2%
GCV = 0.73885 Scale est. = 0.67365 n = 60
#summary(FSG2)
#summary(FSG3)
#summary(FSG4)
summary(FSG1)$r.sq - summary(FSG2)$r.sq
[1] 0.07590024
summary(FSG2)$r.sq - summary(FSG3)$r.sq
[1] 0.05267319
summary(FSG3)$r.sq
[1] 0.0542546
But there is between projects:
ProjGam <- gam(tot_nparticles~ s(depth, k = 3) + factor(project), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500))
summary(ProjGam)
Family: gaussian
Link function: identity
Formula:
tot_nparticles ~ s(depth, k = 3) + factor(project)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.9983 0.4445 20.24 <2e-16 ***
factor(project)P16 17.0001 1.3754 12.36 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1 1 3.552 0.064 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.7 Deviance explained = 70.9%
GCV = 12.411 Scale est. = 11.855 n = 67
library(scales)
Attaching package: ‘scales’
The following object is masked from ‘package:purrr’:
discard
The following object is masked from ‘package:readr’:
col_factor
#https://stackoverflow.com/questions/30179442/plotting-minor-breaks-on-a-log-scale-with-ggplot
log_breaks = function(maj, radix=10) {
function(x) {
minx = floor(min(logb(x,radix), na.rm=T)) - 1
maxx = ceiling(max(logb(x,radix), na.rm=T)) + 1
n_major = maxx - minx + 1
major_breaks = seq(minx, maxx, by=1)
if (maj) {
breaks = major_breaks
} else {
steps = logb(1:(radix-1),radix)
breaks = rep(steps, times=n_major) +
rep(major_breaks, each=radix-1)
}
radix^breaks
}
}
scale_x_log_eng = function(..., radix=10) {
scale_x_continuous(...,
trans=log_trans(radix),
breaks=log_breaks(TRUE, radix),
minor_breaks=log_breaks(FALSE, radix))
}
#theme_set(theme_bw)
PlotPSDmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = psd, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse(limits = c(1200, 0)) + scale_shape_manual(values = c(21:25)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
labs(y = "Depth (m)", x = "Particle Size Distribution Slope") +
geom_hline(yintercept = 175, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue")
#theme_set(theme_cowplot)
PlotParticlesmany <- uds %>%
filter(project == "ETNP") %>%
ggplot(aes(x = tot_nparticles, y = depth, shape = factor(day(time)), fill = hour(time))) +
#geom_path(aes(x = psd_gam)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
geom_point(alpha = .6, size = 2, stroke = 1) +
scale_y_reverse(limits = c(1200, 0)) + scale_shape_manual(values = c(21:25)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
scale_x_log10(breaks = c(10, 100, 1000), minor = c(5, 50, 500)) +
#theme(legend.position = "none") +
#scale_x_log_eng()+
labs(y = "Depth (m)", x = "Particles / L") +
geom_hline(yintercept = 175, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue")
# PlotFluxmany <- uds %>%
# filter(project == "ETNP") %>%
# ggplot(aes(x = tot_flux_fit, y = depth, shape = factor(day(time)), fill = hour(time))) +
#
# #geom_path(aes(x = psd_gam)) +
# #geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1, outline_type = "lower") +
# geom_point(alpha = .6, size = 2, stroke = 1) +
# scale_y_reverse() + scale_shape_manual(values = c(21:25)) +
# scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
# scale_x_log10() + theme(legend.position = "none")
plot_grid(
PlotParticlesmany,
PlotPSDmany,
rel_widths = c(2, 3)
)
Removed 266 rows containing missing values (geom_point).Removed 266 rows containing missing values (geom_point).
ggsave("figures/ParticlesPSDMany.png")
Saving 12 x 7.41 in image
ggsave("figures/ParticlesPSDMany.svg")
Saving 12 x 7.41 in image
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(psd~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG2 <- gam(psd ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG3 <- gam(psd ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG4 <- gam(psd~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.96083 0.01988 -199.3 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.713 1.917 56.197 5.5e-15 ***
s(Day) 1.000 1.000 1.078 0.3036
s(Hour) 1.604 2.000 3.049 0.0285 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.635 Deviance explained = 66.1%
GCV = 0.026006 Scale est. = 0.023702 n = 60
#summary(FSG2)
#summary(FSG3)
summary(FSG4)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.96083 0.01995 -198.5 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.715 1.919 55.141 7.66e-15 ***
s(Hour) 1.565 2.000 2.497 0.0488 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.632 Deviance explained = 65.2%
GCV = 0.02572 Scale est. = 0.023885 n = 60
summary(FSG1)$r.sq - summary(FSG2)$r.sq
[1] 0.03804023
summary(FSG2)$r.sq - summary(FSG3)$r.sq
[1] -0.004905911
summary(FSG3)$r.sq
[1] 0.6015425
Not a significant difference in PSD with respect to time.
But there is between projects:
ProjGam <- gam(psd~ s(depth, k = 3) + factor(project), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500))
summary(ProjGam)
Family: gaussian
Link function: identity
Formula:
psd ~ s(depth, k = 3) + factor(project)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.96167 0.02491 -159.029 <2e-16 ***
factor(project)P16 -0.19977 0.07708 -2.592 0.0118 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.224 1.398 34.54 2.26e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.437 Deviance explained = 45.6%
GCV = 0.039117 Scale est. = 0.037234 n = 67
I wonder if I can show that the profiles aren’t statistically significanlty different. Or that they are for that matter… I think in that case, I run a gam with and without a parameter for profile… And then quantify the effect size of that parameter
Or follow this Gavin Simpson Post https://fromthebottomoftheheap.net/2017/10/10/difference-splines-i/
or anova.gam {mgcv}
Calculate gams for each profile, and then run anova.gam to see if they are different…
PlotNParticlesEP <- uds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(x = tot_nparticles, y = depth, col = project, shape = project)) +
geom_point(alpha = 0.7, size = 2, stroke = 1) +
#geom_path(aes(x = tot_nparticles)) +
#geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + scale_color_manual(values = c("gray20", "brown")) +
labs(x = "Particles/L", y = "Depth (m)") +
theme(legend.position = "none") +
scale_shape_manual(values = c(1:5)) +
geom_hline(yintercept = 175, color = "darkgreen") +
geom_hline(yintercept = 200, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue")
PlotNParticlesEP
I removed one outlyer from p16 for visualization purposes (300 particles/l at surface)
PlotPSDEP <- uds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(x = psd, y = depth, col = project, shape = project)) +
geom_point(alpha = 0.7, size = 2, stroke = 1) +
geom_path(aes(x = psd_gam)) +
geom_ribbon(aes(x = psd_gam, xmin = psd_gam - 2 * psd_seg, xmax = psd_gam + 2 * psd_seg), alpha = 0.1) +
scale_y_reverse(limits = c(1000, 0)) + scale_color_manual(values = c("gray20", "brown")) +
scale_shape_manual(values = c(1:5)) + labs(y = "", x = "Particle Size Distribution Slope") +
geom_hline(yintercept = 175, color = "darkgreen") +
geom_hline(yintercept = 200, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue")
PlotPSDEP
I may just cow these togther.
plot_grid(PlotNParticlesEP, PlotPSDEP, rel_widths = c(2,3), labels = c("A", "B"))
Removed 1211 rows containing missing values (geom_point).Removed 1211 rows containing missing values (geom_point).Removed 1211 row(s) containing missing values (geom_path).
ggsave("figures/ParticlesAndPSD_ETNPVsP16.svg")
Saving 10 x 4 in image
ggsave("figures/ParticlesAndPSD_ETNPVsP16.png")
Saving 10 x 4 in image
mainParticleComponents <- bds %>%
filter(profile %in% c("stn_043", "p16n_100")) %>%
select(project, profile, depth,
tot_nparticles, small_nparticles, big_nparticles,
tot_psd = psd, small_psd, big_psd,
tot_flux_fit, small_flux_fit, big_flux_fit) %>%
pivot_longer(cols = -c("project", "profile", "depth")) %>%
separate(name, c("size", "meas")) %>%
mutate(meas = recode(meas, nparticles = "particles/L")) %>%
mutate(meas = factor(meas, levels = c("particles/L", "flux", "psd")))
Expected 2 pieces. Additional pieces discarded in 273 rows [7, 8, 9, 16, 17, 18, 25, 26, 27, 34, 35, 36, 43, 44, 45, 52, 53, 54, 61, 62, ...].
PlotFlx <- mainParticleComponents %>%
filter(meas != "psd") %>%
ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size ~ meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) + scale_x_log10() + theme(axis.title.x = element_blank(), legend.position = "none", strip.background.y = element_blank(), strip.text.y = element_blank(), plot.margin = unit(c(7,0,7,7), "pt")) + scale_color_manual(values = c("brown", "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = 175, color = "darkgreen")
PlotPSD <- mainParticleComponents %>%
filter(meas == "psd") %>%
ggplot(aes(y = depth, x = value, col = project, shape = project)) + facet_grid(size~meas, scales = "free_x") + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) +
theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.line.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank(), plot.margin = unit(c(7,7,26.5,0), "pt")) +
scale_color_manual(values = c("brown", "gray20")) + scale_shape_manual(values = c(1:5)) + theme(axis.text.x = element_text(angle = 90)) + geom_hline(yintercept = 175, color = "darkgreen")
plot_grid(PlotFlx, PlotPSD, rel_widths = c(3, 2))
Removed 246 rows containing missing values (geom_point).Removed 123 rows containing missing values (geom_point).
ggsave("figures/BigVsSmall.svg")
Saving 7.29 x 4.5 in image
ggsave("figures/BigVsSmall.png")
Saving 7.29 x 4.5 in image
Flux small and flux tot track so closely because ag > psd. since the size distribution of the flux sould be PSD + ag (psd is negative in this case). Yo ucan see the variance at the one depth where psd is flatest at the very top. time 2017-01-13 17:51:31 202.4 L
eg_dataline <- bds %>%
filter(profile == "stn_043", depth == 162.5)
eg_slope = eg_dataline %>% pull(psd)
eg_icp = eg_dataline %>% pull(icp)
eg_vol = eg_dataline %>% pull(vol)
eg_datablock <- bes %>%
filter(profile == "stn_043", depth == 162.5)
eg_lb = eg_datablock$lb
eg_binsize = eg_datablock$binsize
eg_nnp = exp(eg_icp + log(eg_lb) * eg_slope)
eg_np = eg_nnp * eg_binsize
eg_tp = eg_np * eg_vol
eg_df <- tibble(lb = eg_lb, n_nparticles = eg_nnp, nparticles = eg_np, TotalParticles = eg_tp)
EgNNP <- eg_datablock %>%
ggplot(aes(x = lb, y = n_nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs(y = "Binsize & Volume Normalized \n Particles (#/L/mm)", x = "Size (mm)")
EgNP <- eg_datablock %>%
ggplot(aes(x = lb, y = nparticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs(y = "Normalized Particles" , x = "Size (mm)")
EgTP <- eg_datablock %>%
ggplot(aes(x = lb, y = TotalParticles)) + geom_point() + scale_x_log10() + scale_y_log10() +
geom_path(data = eg_df) + labs( y = "Total Particles Observed (#)", x = "Size (mm)")
plot_grid(EgNNP, EgTP, labels = c("A", "B"))
Transformation introduced infinite values in continuous y-axisTransformation introduced infinite values in continuous y-axis
ggsave("figures/ExamplePSD163m.png")
Saving 7.29 x 4.5 in image
ggsave("figures/ExamplePSD163m.svg")
Saving 7.29 x 4.5 in image
bds %>%
ggplot(aes(y = depth, x = Flux_Smooth, col = factor(time))) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(1000, 0)) + scale_x_log10()
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
FSG1 <- gam(Flux_Smooth~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG2 <- gam(Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
FSG3 <- gam(Flux_Smooth ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(FSG1)
Family: gaussian
Link function: identity
Formula:
Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4,
bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.7319 0.9334 37.21 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.875 1.984 11.823 0.000126 ***
s(Day) 1.756 1.939 2.168 0.094012 .
s(Hour) 1.114 2.000 1.658 0.059271 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.33 Deviance explained = 38.4%
GCV = 57.811 Scale est. = 52.277 n = 60
summary(FSG2)
Family: gaussian
Link function: identity
Formula:
Flux_Smooth ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.7319 0.9607 36.15 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.869 1.983 11.102 0.000206 ***
s(Day) 1.811 1.964 2.743 0.094614 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.29 Deviance explained = 33.4%
GCV = 60.062 Scale est. = 55.377 n = 60
summary(FSG3)
Family: gaussian
Link function: identity
Formula:
Flux_Smooth ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 34.7319 0.9944 34.93 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.857 1.979 10.73 0.000277 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.239 Deviance explained = 26.3%
GCV = 62.291 Scale est. = 59.325 n = 60
summary(FSG1)$r.sq - summary(FSG2)$r.sq
[1] 0.03975922
summary(FSG2)$r.sq - summary(FSG3)$r.sq
[1] 0.05063133
summary(FSG3)$r.sq
[1] 0.2392424
bds %>% filter(project == "ETNP") %>% select(profile, depth, Flux_Smooth) %>% pivot_wider(names_from = profile, values_from = Flux_Smooth)
Something is off. All of the flux profiles are identical. Skip this
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
plt1 <- bds %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = DFP, col = factor(time), shape = factor(time))) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(1000, 0)) + xlim(c(0.5, 1.5))+ geom_vline(xintercept = 1) +
scale_color_manual(values = c(rep("black", 5), rep("blue", 5))) + scale_shape_manual(values = rep(1:5, 2))
plotly::ggplotly(plt1)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
What the heck is going on with DFP here. Why is it usually > 1 shouldn’t it be less than 1 when flux is decreasing? This very deep increasing flux seems improbable to me. Lets check the smooths. Or only go to 1000m.
, legend.background = element_blank(), legend.box.background = element_rect()
scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
scientific_10_b <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
scientific_10_c <- function(x) {
xout <- gsub("1e", "10^{", format(x),fixed=TRUE)
xout <- gsub("{-0", "{-", xout,fixed=TRUE)
xout <- gsub("{+", "{", xout,fixed=TRUE)
xout <- gsub("{0", "{", xout,fixed=TRUE)
xout <- paste(xout,"}",sep="")
return(parse(text=xout))
}
scale_x_log10nice <- function(name=NULL,omag=seq(-10,20),...) {
breaks10 <- 10^omag
scale_x_log10(breaks=breaks10,labels=scientific_10_c(breaks10),...)
}
#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlx <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_log10nice()+
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
#labs(x = "moo", y = "Depth (m)") +
geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 20, xmax = 180, ymin = 75, ymax = 500), colour = "red", fill = NA, inherit.aes = FALSE) +
theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm")) +
geom_segment(aes(y = 160, yend = 160, x = 20, xend = 500), color = "darkgreen", stroke = 0.5)+
geom_segment(aes(y = 850, yend = 850, x = 20, xend = 500), color = "darkblue", stroke = 0.5)#+ geom_hline(yintercept = 850, color = "darkblue")
Ignoring unknown parameters: strokeIgnoring unknown parameters: stroke
pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
pltFlxLegend <- get_legend(pltFlx)
Removed 14 rows containing missing values (geom_point).
pltFlx
#plotly::ggplotly(plt1)
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxZoom <- bds %>% filter(project == "ETNP" & depth <= 500 & depth >= 75) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse()+
#scale_x_log10() +
scale_x_log10(breaks = c(seq(from = 20, to = 50, by = 10), seq(from = 60, to = 180, by = 20)), limits = c(20, 180)) +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(values = rep(21:25, 2)) +
scale_fill_gradientn(breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
theme(axis.text.x = element_text(angle = 90)) +
labs(x = "Smoothed Flux", y = "Depth") + theme(legend.position = "none")+
geom_hline(yintercept = 160, color = "darkgreen")
pltFlxZoom
#plotly::ggplotly(plt1)
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3 <- bds %>% filter(project == "ETNP") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2, stroke = 1)+
#geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_continuous(limits = c(-2.1, .6), breaks = seq(from = -2, to = .75, by = 0.5)) +
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
geom_vline(xintercept = 0) +
labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")+
geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = 850, color = "darkblue")
#labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")
pltDelta3
#plotly::ggplotly(plt1pos)
# #plot_grid(pltFlxNoLegend, pltFlxZoom, pltDelta3, pltFlxLegend)
#
# pltFlxLegend <- get_legend(pltFlx + theme(legend.box.margin = margin(0, 0, 40, 10)))
#
# pgTop <- plot_grid(pltFlxNoLegend, pltFlxZoom + theme(plot.margin = unit(c(1, 0, 3, 0), units = "cm")), rel_widths = c(2, 1), labels = c("A", "B"))
# pgBottom <- plot_grid(pltDelta3, pltFlxLegend , rel_widths = c(3, 1), labels = c("C", ""), label_size = 14)
# pgBoth <- plot_grid(pgTop, pgBottom, ncol = 1)
#
# pgBoth
#
# ggsave("figures/FluxDeepDive.png")
# ggsave("figures/FluxDeepDive.svg")
Within panel drawing
pgTop <- ggdraw(pltFlxNoLegend
) +
draw_plot(pltFlxZoom, .4, .25, .55, .60) +
draw_plot_label(
c("","B"),
c(.05, 0.55),
c(1, 0.85),
size = 16
)
Removed 14 rows containing missing values (geom_point).
pgTop
pgBottom <- plot_grid(pltDelta3, pltFlxLegend , rel_widths = c(3, 1), labels = c(“C”, ""), label_size = 14)
I don’t know whats going on below here
pgBottom <- pltDelta3 + geom_rect(aes(xmin = -2, xmax = -1.15, ymin = 170, ymax = 1000), colour = "gray50", fill = NA, inherit.aes = FALSE) + draw_plot(pltFlxLegend , -1.9, -575, .7)
pgBoth <- plot_grid(pgTop + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
pgBottom + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
ncol = 1, rel_heights = c(4, 4), labels = c("A", "C"), label_size = 16)
Removed 33 rows containing missing values (geom_point).
pgBoth
ggsave("figures/FluxDeepDive.png")
Saving 5 x 9 in image
ggsave("figures/FluxDeepDive.svg")
Saving 5 x 9 in image
# #plot_grid(pltFlxNoLegend, pltFlxZoom, pltDelta3, pltFlxLegend)
#
# pltFlxLegend <- get_legend(pltFlx + theme(legend.box.margin = margin(0, 0, 40, 10)))
#
# pgTop <- plot_grid(pltFlxNoLegend + ylim(c(1000, 0)), pltFlxZoom + theme(plot.margin = unit(c(1, 0, 3, 0), units = "cm")), rel_widths = c(2, 1), labels = c("A", "B"))
# pgBottom <- plot_grid(pltDelta3 + ylim(c(1000, 0)), pltFlxLegend , rel_widths = c(3, 1), labels = c("C", ""))
# pgBoth <- plot_grid(pgTop, pgBottom, ncol = 1)
#
# pgBoth
#
# #ggsave("figures/FluxShallowDive.png")
# #ggsave("figures/FluxShallowDive.svg")
Test for day to day and hourly variability in delta flux
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
DFG1 <- gam(pracma::nthroot(DF/DZ, 5)~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG2 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
DFG3 <- gam(pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 250 & depth <=500 & project == "ETNP"))
summary(DFG1)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3) +
s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.05564 -3.305 0.00217 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.722 1.922 4.549 0.03760 *
s(Day) 1.943 1.997 8.215 0.00139 **
s(Hour) 1.470 2.000 2.587 0.04015 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.359 Deviance explained = 43.9%
GCV = 0.15225 Scale est. = 0.13001 n = 42
summary(DFG2)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.05943 -3.094 0.00373 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.661 1.885 3.673 0.07425 .
s(Day) 1.910 1.992 5.174 0.00943 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.268 Deviance explained = 33.2%
GCV = 0.16645 Scale est. = 0.14833 n = 42
summary(DFG3)
Family: gaussian
Link function: identity
Formula:
pracma::nthroot(DF/DZ, 5) ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.18385 0.06635 -2.771 0.00849 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.592 1.834 3.213 0.107
R-sq.(adj) = 0.0878 Deviance explained = 12.3%
GCV = 0.19707 Scale est. = 0.18491 n = 42
summary(DFG1)$r.sq - summary(DFG2)$r.sq
[1] 0.09038471
summary(DFG2)$r.sq - summary(DFG3)$r.sq
[1] 0.1804172
summary(DFG3)$r.sq
[1] 0.08782536
png(filename = “./figures/CombinedP2Info.png”, width = 10, height = 8, units = “in”, res = 200) StationInfoPlot() dev.off()
#plot.new()
FluxGamPlot <- function(){
par(mfrow = c(2,2))
plot(DFG1)
mtext(expression(bold("C")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,1))
mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,2))
mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
}
FluxGamPlot()
png(filename = "./figures/FluxGamPlot.png", width = 10, height = 8, units = "in", res = 200)
FluxGamPlot()
dev.off()
png
2
Check of actual data for hour
ggplot(data = bds %>% filter(depth >= 175, depth <= 500), aes(y = DF/DZ, x = hour(time), col = depth, group = depth)) + geom_point() + geom_line()
#Osps
(u mol C / m^3 / day)
disagFig <- bds %>% filter(project == "ETNP") %>%
ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), shape = factor(day(time)), fill = hour(time), group = factor(time))) + geom_point(size = 2) + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(-1, 1)) +
geom_vline(xintercept = 0) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
#labs(x = bquote("Observed - Modeled Small Particle Flux"~(μmol/m^3/day)), y = "Depth (m)") +
labs(x = paste("Deviation from Model", expression((μmol/m^3/day)))) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 850, color = "darkblue")
disagFig
#ggsave("..figures/FluxSizeShift.svg"
ggsave("figures/FluxSizeShift.png")
Saving 6 x 4 in image
ggsave("figures/FluxSizeShift.svg")
Saving 6 x 4 in image
For NSF Proposal plot_grid( PlotParticlesmany, PlotPSDmany, rel_widths = c(2, 3) )
pgTop <- ggdraw(pltFlxNoLegend
) +
draw_plot(pltFlxZoom, .4, .25, .55, .60) +
draw_plot_label(
c("","D"),
c(.05, 0.55),
c(1, 0.85),
size = 16
)
Removed 14 rows containing missing values (geom_point).
pgCF <- plot_grid(PlotParticlesmany + ylim(c(1000, 0)),
PlotPSDmany + ylim(c(1000, 0)) +theme(legend.position = "none"),
pgTop + theme(plot.margin = unit(c(0, 0, 0, 0), units = "cm")),
disagFig + theme(plot.margin = unit(c(0, 0, .2, 0), units = "cm")),
ncol = 2, rel_heights = c(4, 4), labels = c("A", "B", "C", "E"), label_size = 16)
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
Removed 308 rows containing missing values (geom_point).Removed 308 rows containing missing values (geom_point).Removed 41 rows containing missing values (geom_point).
pgCF
ggsave("./figures/UVP_for_zoop_proposal.png", pgCF, width = 10, height = 8, units = "in")
bdsAddTime <- bds %>%
mutate(Hour = hour(time), Day = day(time))
OZG1 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc"), knots = list(Hour = c(0, 24)), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
OZG2 <- gam(ospsDZ ~ s(depth, k = 3) + s(Day, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
OZG3 <- gam(ospsDZ ~ s(depth, k = 3), data = bdsAddTime %>% filter(depth >= 175 & depth <=500 & project == "ETNP"))
summary(OZG1)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3) + s(Hour, k = 4, bs = "cc")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.066685 0.009465 7.046 3.43e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.681 1.898 16.665 1.66e-05 ***
s(Day) 1.870 1.982 7.601 0.000825 ***
s(Hour) 1.319 2.000 2.124 0.047647 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.45 Deviance explained = 49.6%
GCV = 0.0059577 Scale est. = 0.0053749 n = 60
summary(OZG2)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3) + s(Day, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.066685 0.009813 6.796 7.85e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.644 1.873 16.338 2.28e-05 ***
s(Day) 1.904 1.991 7.155 0.00258 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.409 Deviance explained = 44.5%
GCV = 0.0062509 Scale est. = 0.0057771 n = 60
summary(OZG3)
Family: gaussian
Link function: identity
Formula:
ospsDZ ~ s(depth, k = 3)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.06669 0.01084 6.154 7.82e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(depth) 1.547 1.794 14.83 7.66e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.279 Deviance explained = 29.8%
GCV = 0.0073573 Scale est. = 0.007045 n = 60
summary(OZG1)$r.sq - summary(OZG2)$r.sq
[1] 0.04114558
summary(OZG2)$r.sq - summary(OZG3)$r.sq
[1] 0.1296981
summary(OZG3)$r.sq
[1] 0.2793194
OSMSGamPlot <- function(){
par(mfrow = c(1,2))
plot(OZG2)
mtext(expression(bold("B")), side = 3, line = 0, adj = 0, cex = 2)
par(mfg = c(1,1))
mtext(expression(bold("A")), side = 3, line = 0, adj = 0, cex = 2)
}
OSMSGamPlot()
png(filename = "./figures/OSMSGamPlot.png", width = 10, height = 6, units = "in", res = 200)
OSMSGamPlot()
dev.off()
png
2
plot(OZG2)
bds %>% filter(project == "P16") %>%
ggplot(aes(y = depth, x = ospsDZ)) + facet_wrap(~project) + geom_point() + scale_y_reverse(limits = c(500, 0)) + geom_vline(xintercept = 0)
trapFlux3 <- read_csv("dataOut/fluxMS_distilled.csv")
Parsed with column specification:
cols(
Class = [31mcol_character()[39m,
Depth = [32mcol_double()[39m,
TrapID = [31mcol_character()[39m,
TrapType = [31mcol_character()[39m,
SampleType = [31mcol_character()[39m,
C_flux = [32mcol_double()[39m,
C_flux_umol = [32mcol_double()[39m
)
UVPFluxComb <- read_csv("dataOut/CombinedProfileFluxEst_DS.csv")
Parsed with column specification:
cols(
depth = [32mcol_double()[39m,
Flux = [32mcol_double()[39m
)
UVPFluxOE <- read_csv("dataOut/ObservedVsExpectedFlux.csv")
Parsed with column specification:
cols(
depth = [32mcol_double()[39m,
tn_flux = [32mcol_double()[39m,
profile = [31mcol_character()[39m,
project = [31mcol_character()[39m,
time = [31mcol_character()[39m,
tot_flux2 = [32mcol_double()[39m
)
trapFlux3
UVPFluxComb
fluxMS_distilled_toPlot <- trapFlux3 %>%
mutate(SampleType = recode(SampleType, `plus.p` = "plus-particles", top = "top-collector"))
Remove traps where mass spec didn’t work correctly 2-17 150 1-12 73m 1-12 148 2-14 100 |(TrapID == “2-17” & Depth == 150)
fluxMS_distilled_toPlot2 <- fluxMS_distilled_toPlot %>%
filter(!((TrapID == "1-12") | (TrapID == "2-14" & Depth == 100)|(TrapID == "2-17" & Depth == 150)))
fluxMS_distilled_toPlot2
# UVPFluxComb %>%
# ggplot(aes(y = depth)) + scale_y_reverse(limits = c(1000, 0)) +
# scale_x_continuous(limits = c(0, 200)) +
# geom_point(aes(y = Depth, x = C_flux_umol, fill = SampleType, shape = TrapType),
# colour = "black", stroke = 1, size = 5, data = fluxMS_distilled_toPlot) +
# geom_point(aes(x = Flux), size = 3, shape = 21, color = "white", fill = "black") +
# geom_point(aes(x = -1, y = -1, size = "UVP")) + # dummy point for the legend
# scale_shape_manual(values = c(25, 22))+
# scale_size_manual(values = 1, name = "") +
# ylab("Depth (m)") + xlab("Flux µmolC/m^2/day") +
# guides(fill = guide_legend(override.aes = list(shape = 21))) +
# scale_fill_viridis_d() +
# theme_cowplot() +
# theme(
# legend.position = c(0.5, 0.4),
# legend.box.background = element_rect(color = "black", size = 0.5),
# legend.margin = margin(-10, 5, 10, 5)
# ) +
# geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 15, xmax = 32, ymin = 45, ymax = 195), colour = "red", fill = NA, inherit.aes = FALSE)
# ggsave("figures/FittedFlux.png")
# ggsave("figures/FittedFlux.svg")
Traps where mass spec didn’t work.
UVPFluxPlot00 <- UVPFluxComb %>%
ggplot(aes(y = depth)) + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(0, 200)) +
geom_point(aes(y = Depth, x = C_flux_umol, shape = TrapType, ID = TrapID),
colour = "black", stroke = 1, size = 5, data = fluxMS_distilled_toPlot2) +
geom_line(aes(x = Flux), size = 1, color = "black") +
geom_point(aes(x = -1, y = -1, size = "UVP Estimate")) + # dummy point for the legend
geom_point(aes(x = tot_flux2), size = 3, shape = 21, color = "white", fill = "black", data = UVPFluxOE) +
scale_shape_manual(values = c(25, 22))+
scale_size_manual(values = 1, name = "") +
ylab("Depth (m)") +
#xlab(expression(Flux µmolC/m^2/day)) +
xlab(expression(paste("x axis ", ring(A)^2))) +
xlab(expression(paste("Flux (µ mol C/", m^2, "/day)"))) +
guides(fill = guide_legend(override.aes = list(shape = 21))) +
theme_cowplot() +
theme(
legend.position = c(0.5, 0.4),
legend.box.background = element_rect(color = "black", size = 0.5),
legend.margin = margin(-10, 5, 10, 5)
)
Ignoring unknown aesthetics: ID
# UVPFluxPlot <- UVPFluxPlot00 +
# geom_rect(data = data.frame(project = "ETNP"), aes(xmin = 15, xmax = 32, ymin = 45, ymax = 195), colour = "red", fill = NA, inherit.aes = FALSE)
UVPFluxPlot00
ggsave("figures/FittedFlux.png")
Saving 7.29 x 4.5 in image
ggsave("figures/FittedFlux.svg")
Saving 7.29 x 4.5 in image
horizontalGamPlot <- dataGamHorizontal %>% ggplot(aes(x = resp_fit, y = depth, col = log(lb), group = lb)) + scale_y_reverse() + geom_point() + scale_x_log10(limits = c(10^-8, NA)) + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + geom_errorbar(aes(xmin = resp_lower, xmax = resp_upper), width = 10, alpha = 0.5)+ theme_bw()
TPPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = TotalParticles, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "TotalParticles Observed (#)")
nnpPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = n_nparticles, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "Binsize and Volume Normalized Particles (#/L/mm)")
FitPlot <- bes %>% filter(profile == "stn_043") %>% group_by(lb) %>% ggplot(aes(x = nnp_smooth, xmin = nnp_lower, xmax = nnp_upper, y = depth, col = log(lb), group = lb)) + scale_y_reverse(limits = c(1000, 0)) + geom_point() + scale_x_log10() + scale_color_viridis_c() + geom_path() + geom_vline(xintercept = 1) + geom_vline(xintercept = 5) + labs(y = "Depth (m)", x = "Smoothed - Normalized Particles (#/L/mm)") + geom_errorbar(width = 10, alpha = 0.5)
npLegend <- get_legend(FitPlot + theme(legend.box.margin = margin(0, 0, 40, 200)) + labs(col = expression(log[e](Size (mm)))))
Removed 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).
plot_grid(
TPPlot + theme(legend.position = "none"),
nnpPlot + theme(legend.position = "none"),
npLegend ,
FitPlot + theme(legend.position = "none")
)
Transformation introduced infinite values in continuous x-axisTransformation introduced infinite values in continuous x-axisRemoved 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).Transformation introduced infinite values in continuous x-axisTransformation introduced infinite values in continuous x-axisRemoved 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).Removed 325 rows containing missing values (geom_point).Removed 325 row(s) containing missing values (geom_path).
ggsave("figures/AllParticleSizes.svg")
Saving 10 x 6.18 in image
ggsave("figures/AllParticleSizes.png")
Saving 10 x 6.18 in image
SameGam <- gam(TotalParticles ~s(log(lb), log(depth), by = factor(time)), offset = log(vol * binsize), family = nb(),
data = bes %>% filter(project == "ETNP"))
besE <- bes %>% filter(project == "ETNP")
lb_new <- exp(seq(from = log(0.1), to = log(2.1), by = 0.05))
ub_new <- lead(lb_new)
binsize_new <- ub_new - lb_new
lbbs <- tibble(lb = lb_new, ub = ub_new, binsize = binsize_new)
Expanded <- expand_grid(lb = exp(seq(from = log(0.1), to = log(2), by = 0.05)), depth = seq(from = 20, to = 2000, by = 20), time = as.factor(unique(besE$time))) %>% left_join(lbbs, by = "lb")
Pred <- exp(predict(SameGam, Expanded))
ToPlot <- bind_cols(Expanded, nnparticles = Pred) %>% mutate(time = as.character(time)) %>% mutate(nparticles = nnparticles * binsize)
ToPlot %>% filter(lb <= 2) %>% ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c() + scale_y_reverse() + scale_x_log10() + facet_wrap(~time) + geom_contour(color = "black")
meanBese <- ToPlot %>% filter(lb <= 2) %>% group_by(lb, depth) %>% summarize(nparticles = mean(nparticles), nnparticles = mean(nnparticles))
WBColorMap <- meanBese%>%
ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c(name = "log10(number density \n (normalized))") + scale_y_reverse() + scale_x_log10() + geom_contour(color = "black") + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = 850, color = "darkblue")
WBColorMap
Average of everything
#meanBese043 <- ToPlot %>% filter(lb <= 2, time == "2017-01-13 11:51:31")
meanBese%>%
ggplot(aes(x = lb, y = depth, fill = log10(nnparticles), z = log10(nnparticles))) + geom_tile() + scale_fill_viridis_c() + scale_y_reverse() + scale_x_log10() + geom_contour(color = "black") + geom_hline(yintercept = 160, color = "darkgreen")
Just 043
mbGam <- meanBese %>% group_by(depth) %>% nest() %>%
mutate(mod = map(data, ~gam(log(nnparticles) ~ log(lb), family = gaussian(), data = .))) %>%
mutate(psd = map_dbl(mod, ~summary(.)$p.coeff[2]))
mbGam %>% ggplot(aes(x = psd, y = depth)) + geom_path() + scale_y_reverse() + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = 850, color = "darkblue")
mbGam <- meanBese043 %>% group_by(depth) %>% nest() %>% mutate(mod = map(data, ~gam(log(nnparticles) ~ log(lb), family = gaussian(), data = .))) %>% mutate(psd = map_dbl(mod, ~summary(.)$p.coeff[2]))
pWBPSD <- mbGam %>% ggplot(aes(x = psd, y = depth)) + geom_path() + scale_y_reverse() + geom_hline(yintercept = 160, color = "darkgreen") + geom_hline(yintercept = 850, color = "darkblue")
pWBPSD
bds %>% filter(profile == “stn_043”, depth <= 2000) %>% ggplot(aes(x = psd_gam, xmin = psd_gam - psd_seg * 2, xmax = psd_gam + psd_seg * 2, y = depth)) + geom_path(size = 1) + scale_y_reverse() + geom_hline(yintercept = 175, color = “darkgreen”) + geom_hline(yintercept = 950, color = “darkblue”) + geom_ribbon(alpha = 0.2) + labs(x = “PSD slope”)
All of them
bds %>% filter(profile == "stn_043", depth <= 2000) %>% ggplot(aes(x = psd_gam, xmin = psd_gam - psd_seg * 2, xmax = psd_gam + psd_seg * 2, y = depth)) + geom_path(size = 1) + scale_y_reverse() + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue") + geom_ribbon(alpha = 0.2) + labs(x = "PSD slope")
043 only
bds %>% filter(profile == "stn_043", depth <= 2000, depth > 175) %>% ggplot(aes(x = small_biovolume, y = depth)) + geom_path(size = 1) + scale_y_reverse() + geom_hline(yintercept = 175, color = "darkgreen") + geom_hline(yintercept = 950, color = "darkblue") + geom_point()
ubDf0 <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global)
ubDf <- ubDf0 %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth)
photicBiomass <- ubDf %>% filter(depth <= 180, depth >= 160) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
ubDf <- ubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length
ubDf %>% ggplot(aes(x = nbiomass, y = depth , group = time, col = time)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1))
ubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- ubDf %>% filter(depth <= 180, depth >= 160) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
ubDf <- ubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
ubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1)) + geom_hline(yintercept = 175, color = "darkgreen")
PubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb < 0.5) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- PubDf %>% filter(depth <= 165, depth >= 155) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
PubDf <- PubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBS <- PubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1.2)) + geom_hline(yintercept = 160, color = "darkgreen") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = 850, color = "darkblue") + labs( x = "Small particle mass (norm.)")
pWBS
LubDf <- ToPlot %>% mutate(ubiomass = nparticles * lb ^ ag_global) %>% filter(lb >= 0.5) %>% group_by(time, depth) %>% summarize(ubiomass = sum(ubiomass)) %>% ungroup %>% group_by(depth) %>% summarise(ubiomass = mean(ubiomass))
photicBiomass <- LubDf %>% filter(depth <= 165, depth >=155) %>% summarize(ubiomass = mean(ubiomass)) %>% pull(ubiomass)
LubDf <- LubDf %>% mutate(nbiomass = ubiomass/photicBiomass)
pWBL <- LubDf %>% ggplot(aes(x = nbiomass, y = depth)) + geom_path() + scale_y_reverse() + scale_x_continuous(limits = c(0,1)) + geom_hline(yintercept = 160, color = "darkgreen") + labs( x = "Large particle mass (norm.)") + geom_vline(xintercept = 1, color = "gray50") + geom_vline(xintercept = 0, color = "gray50") + geom_hline(yintercept = 850, color = "darkblue")
pWBL
For tom and danielle
WBColorMap
pWBPSD
pWBS
pWBL
WBFig5 <- plot_grid(pWBPSD, pWBS,pWBL, nrow = 1, labels = c("B", "C", "D"))
Removed 6 row(s) containing missing values (geom_path).Removed 7 row(s) containing missing values (geom_path).
WBFig5
WBcombined <- plot_grid(WBColorMap + theme(plot.margin = unit(c(0,3,0, 3), "cm")), WBFig5, ncol = 1, labels = c("A", ""))
WBcombined
ggsave("figures/WBModelValidation.png")
Saving 8 x 6 in image
scientific_10 <- function(x) {parse(text=gsub("e\\+*", " %*% 10^", scales::scientific_format()(x))) }
#https://stackoverflow.com/questions/10762287/how-can-i-format-axis-labels-with-exponents-with-ggplot2-and-scales
#jacob_magnitude <- function(x){expression(10^round(log10(x)))}
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltFlxP16 <- bds %>% filter(project == "P16") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = Flux_Smooth, group = factor(time))) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_log10(limits = c(35, 150),breaks = seq(from = 20, to = 150, by = 20)) +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
labs(x = bquote(Smoothed~Flux~(µmol~C/m^2/d)), y = "Depth (m)") +
geom_hline(yintercept = 200, color = "darkgreen") +
theme(axis.text.x = element_text(angle = 90, vjust = .3), legend.spacing = unit(.1, "cm"))
#
#
#
# pltFlxNoLegend <- pltFlx + theme(legend.position = "none")
# pltFlxLegend <- get_legend(pltFlx)
#
pltFlxP16
# #plotly::ggplotly(plt1)
cb10 <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
pltDelta3P16 <- bds %>% filter(project == "P16") %>% #filter(DFP > 1) %>% #filter(profile %in% c("stn_043", "p16n_100")) %>%
ggplot(aes(y = depth, x = pracma::nthroot(DF/DZ, 5), group = factor(time))) + geom_point(size = 3, stroke = 1)+
geom_path() +
scale_y_reverse(limits = c(1000, 0))+
scale_x_continuous(limits = c(-1, .1), breaks = seq(from = -2, to = .75, by = 0.5)) +
#scale_x_log10() +
scale_color_gradient2(low = "darkgreen", mid = "gray80", high = "purple", midpoint = 10) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) +
scale_fill_gradientn(name = "Hour", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) +
geom_vline(xintercept = 0) +
geom_hline(yintercept = 200, color = "darkgreen")+
labs(x = bquote((DF/DZ)^{1/5}~(µmolC/m^3/d)^{1/5}), y = "Depth (m)") + theme(legend.pos = "none")
#labs(x = "(DF/DZ) ^ 1/5 (µmol C/m^3/d) ^ 1/5")
pltDelta3P16
#plotly::ggplotly(plt1pos)
osms_p16 <- bds %>% filter(project == "P16") %>%
ggplot(aes(y = depth, x = pracma::nthroot(ospsDZ, 3), group = factor(time))) + geom_point(size = 3) + geom_path() + scale_y_reverse(limits = c(1000, 0)) +
scale_x_continuous(limits = c(-1, 1)) +
geom_vline(xintercept = 0) + scale_shape_manual(name = "Day of Month", values = rep(21:25, 2)) + labs(x = "Observed - Modeled Small Particle Flux \n µmol/m^3/day") +
scale_fill_gradientn(name = "Hour of Day", breaks = c(0, 6, 12, 18, 24), colors = c("black", "blue", "white", "orange", "black")) + geom_hline(yintercept = 175, color = "darkgreen")
plotly::ggplotly(osms_p16)
#ggsave("..figures/FluxSizeShift.svg"
plot_grid(
pltFlxP16,
pltDelta3P16,
osms_p16
)
Removed 28 rows containing missing values (geom_point).Removed 28 row(s) containing missing values (geom_path).Removed 32 rows containing missing values (geom_point).Removed 32 row(s) containing missing values (geom_path).Removed 29 rows containing missing values (geom_point).Removed 29 row(s) containing missing values (geom_path).
ggsave("figures/P16FluxRelate.svg")
Saving 8 x 8 in image
ggsave("figures/P16FluxRelate.png")
Saving 8 x 8 in image
source("ModelStuff.R")
scan_for_example <- bds %>% filter(project == "ETNP", depth < 500, depth > 200) %>% select(profile, depth, DFP, use_this_DFP, ospsDZ)
#loc_station = "stn_036"
loc_station = "stn_043"
loc_depth = 225
loc_prev_depth = 112.5
allDFPs <- bds %>% filter(profile == loc_station, depth >= loc_prev_depth, depth <= loc_depth) %>% summarize(DFP = prod(DFP), use_this_DFP = prod(use_this_DFP))
loc_DFP <- allDFPs %>% pull(DFP)
loc_use_DFP <- allDFPs %>% pull(use_this_DFP)
for_single_disag <- bes %>% filter(profile == loc_station, depth %in% c(loc_prev_depth, loc_depth)) %>% select(depth, lb, nnp_smooth) %>%
mutate(depth = recode(depth, `112.5` = "Shallow", `225` = "Deep")) %>% # I have no idea how to not hard code this bit
pivot_wider(names_from = depth, values_from = nnp_smooth)
with_disag <- for_single_disag %>%
mutate(Predicted_Deep = remin_smooth_shuffle(Shallow, loc_use_DFP))
#remin_smooth_shuffle(for_single_disag$Shallow,loc_use_DFP)
for_plot_disag <- with_disag %>% pivot_longer(cols = -lb) %>% #filter(lb <= 5) %>%
mutate(name = factor(name, levels = c("Shallow", "Deep", "Predicted_Deep"))) %>%
mutate(name = recode_factor(name, Shallow = "Shallow (112.5m)", Deep = "Deep (225m)", Predicted_Deep = "Predicted Deep (225m)"))
for_plot_disag %>% ggplot(aes(x = lb, y = value, shape = name)) + geom_point() + scale_x_log10() + scale_y_log10() + scale_shape_manual(values = c(1, 6, 3)) + theme(legend.title = element_blank()) + labs(x = "Particle Size (mm)", y = "Normalized Particle Abundance (#/L/mm)")
ggsave("figures/DisagExample.png")
Saving 7.29 x 4.5 in image
ggsave("figures/DisagExample.svg")
Saving 7.29 x 4.5 in image
dataBinned <- read_csv("data/backscatter_table_go7.csv")
Parsed with column specification:
cols(
frequency = [32mcol_double()[39m,
depth_bin = [32mcol_double()[39m,
time_bin = [34mcol_datetime(format = "")[39m,
value = [32mcol_double()[39m
)
dataBinned_01 <- dataBinned %>%
mutate(timeMex = with_tz(time_bin, tzone = "US/Central") )
startDay <- dataBinned_01$timeMex %>% na.omit %>% min %>% floor_date(unit = "days")
endDay <- dataBinned_01$timeMex %>% na.omit %>% max %>% ceiling_date(unit = "days")
timeBreaks <- seq(from = startDay, to = endDay, by = "12 hours")
timeLabels <- format(timeBreaks)
plotLetters <- tribble(
~letter, ~depth_bin, ~timeMex,
"A", 350, as.POSIXct("2017-01-07 13:00:00"),
"B", 200, as.POSIXct("2017-01-09 23:00:00"),
"C", 150, as.POSIXct("2017-01-08 13:00:00"),
"D", 625, as.POSIXct("2017-01-12 07:00:00"),
"E", 750, as.POSIXct("2017-01-10 02:00:00")
)
library(shadowtext)
plot18k <- dataBinned_01 %>% filter(frequency == 18000) %>% ggplot(aes(x = timeMex, y = depth_bin, fill = value)) + geom_tile() + scale_y_reverse() + scale_fill_viridis_c(limits = c(-165, -75), oob = scales::squish) +
scale_x_datetime(breaks = timeBreaks, date_labels = "%d::%H") + labs(x = "day::hour", y = "depth (m)", fill = "backscatter (dB)") + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + geom_hline(yintercept = 160, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue") +
geom_shadowtext(data = plotLetters, aes(x = timeMex - 12 * 60^2, y = depth_bin, label = letter), inherit.aes = FALSE, size = 6, bg.color = "white", color = "black") +
geom_segment(data = plotLetters, inherit.aes = FALSE, aes(x = timeMex - 9 * 60^2, xend = timeMex, y = depth_bin, yend = depth_bin), arrow = arrow(length = unit(0.03, "npc")), color = "white", size = 1.5) +
geom_segment(data = plotLetters, inherit.aes = FALSE, aes(x = timeMex - 9 * 60^2, xend = timeMex, y = depth_bin, yend = depth_bin), arrow = arrow(length = unit(0.03, "npc")))
plot18k
#ggsave("figures/stationP2_EK60_18kOnly.png")
plotLetters <- tribble(
~letter, ~depth_bin, ~timeMex,
"A", 300, as.POSIXct("2017-01-07 13:00:00"),
"B", 200, as.POSIXct("2017-01-09 23:00:00"),
"C", 150, as.POSIXct("2017-01-08 13:00:00"),
"D", 625, as.POSIXct("2017-01-12 07:00:00"),
"E", 750, as.POSIXct("2017-01-10 02:00:00")
)
library(shadowtext)
plot38k <- dataBinned_01 %>% filter(frequency == 38000) %>% ggplot(aes(x = timeMex, y = depth_bin, fill = value)) + geom_tile() + scale_y_reverse() + scale_fill_viridis_c(limits = c(-165, -95), oob = scales::squish) +
scale_x_datetime(breaks = timeBreaks, date_labels = "%d::%H") + labs(x = "Day::Hour", y = "Depth (m)", fill = "Backscatter (dB)") + theme_bw() + geom_hline(yintercept = 160, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue") +
geom_shadowtext(data = plotLetters[1:2,], aes(x = timeMex - 12 * 60^2, y = depth_bin, label = letter), inherit.aes = FALSE, size = 6, bg.color = "white", color = "black") +
geom_segment(data = plotLetters[1:2,], inherit.aes = FALSE, aes(x = timeMex - 9 * 60^2, xend = timeMex, y = depth_bin, yend = depth_bin), arrow = arrow(length = unit(0.03, "npc")), color = "white", size = 1.5) +
geom_segment(data = plotLetters[1:2,], inherit.aes = FALSE, aes(x = timeMex - 9 * 60^2, xend = timeMex, y = depth_bin, yend = depth_bin), arrow = arrow(length = unit(0.03, "npc"))) +
theme(axis.text.x = element_text(size = 12,angle = 90, vjust = 0.5),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 14),
legend.title = element_text(size = 14)
)
plot38k
ggsave("figures/stationP2_EK60_38kOnly.png")
Saving 7.29 x 4.5 in image
plot18k + scale_y_reverse(limits = c(500, 200))
Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.
plot200k <- dataBinned_01 %>% filter(frequency == 200000) %>% ggplot(aes(x = timeMex, y = depth_bin, fill = value)) + geom_tile() + scale_y_reverse() + scale_fill_viridis_c(limits = c(-165, -155), oob = scales::squish) +
scale_x_datetime(breaks = timeBreaks, date_labels = "%d::%H") + labs(x = "day::hour", y = "depth (m)", fill = "backscatter (dB)") + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) + geom_hline(yintercept = 160, color = "darkgreen") +
geom_hline(yintercept = 850, color = "darkblue")
plot200k
dataBinned_01 %>% ggplot(aes(x = timeMex, y = depth_bin, fill = value)) + geom_tile() + scale_y_reverse() + scale_fill_viridis_c(limits = c(-165, -75), oob = scales::squish) +
scale_x_datetime(breaks = timeBreaks, date_labels = "%d::%H") + labs(x = "day::hour", y = "Depth (m)", fill = "backscatter (dB)") + theme_bw() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
facet_wrap(~frequency, ncol = 1)
#ggsave("figures/stationP2_EK60_go7.svg", width = 4, height = 10)
ggsave("figures/stationP2_EK60_go7.png", width = 6, height = 10)